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1.  INTRODUCTION 


Factors  which  contribute  to  the  dissipation  of  acoustic  energy  in 
the  chamber  of  solid  rocket  motors  during  combustion  include  -  nozzle 
performance,  particulate  matter,  propellant  structural  response, 
rotational  flow,  vortex  shedding,  etc.  The  most  significant  effect 
appears  to  be  the  response  of  the  propellant  combustion  zone  to  acoustic 
pressure  and  acoustic  velocity  oscillations.  A  large  body  of  literature 
exists  relative  to  this  subject,  the  study  of  which  has  been  pioneered 
by  Crocco  [1],  Cantrell  and  Hart  [2],  Culick  (3,4],  and  others.  Flandro 
and  Jacobs  [5],  among  others,  have  noted  that  vortex  shedding  can  also 
lead  to  an  excitation  of  pressure  oscillations  in  solid  propellant 
rocket  motors.  It  is  also  quite  possible  that  high  speed  mean  flows 
affect  the  stability  [6]  significantly.  Effects  of  transonic  flow, 
shock  waves,  fluid  viscosity  at  shear  layers,  turbulence,  nonlinearity 
(second  order  perturb^ -i  ons)  ,  and  radiation  through  high  temperature 
should  also  be  of  concern. 

The  topics  of  study  in  this  paper  are  limited  to  the  basic 
question  of  correct  stability  integral  formulation,  and  to  three- 
dimensional  finite  element  applications  for  stability  calcuations  built 
upon  the  earlier  work  of  Hackett  [7] .  Because  of  the  flexibility  of  the 
finite  element  formulations,  the  present  work  can  be  extended,  without 
difficulty,  to  a  more  general  case  incorporating  various  effects  such 
as  high  speed  flow,  shock  waves,  oarticle  and  structural  damping,  turbu¬ 
lence,  and  radiation. 

In  what  follows  we  present  the  formulation  which  includes  con¬ 
sideration  of  vortex  shedding  and  fluid  viscosity.  It  (has  been  shown 
that  a  simple  and  special  case  of  the  general  formulation  reduces  to 
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the  well-known  results  such  as  "flow  turning"  as  well  as  "velocity 
and  pressure  coupling"  terms.  Numerical  results  of  some  simple  geometries 
are  discussed,  pending  a  full  scale  computer  code  development  based  on  the 

t 

present  formulation. 

2.  FORMULATION  AND  GOVERNING  EQUATIONS 

Consider  a  compressible,  viscous  fluid  and  the  corresponding  non- 
dimensional  equations  for  continuity,  momentum,  energy,  and  state,  as 
follows : 


f?  +  (pui),i 


^ui  1  1  /  1  \ 

P  —  +  pu .  .  U  .4 - p  .  --  rr—  fu.  ..+■=•  U.  .  .  )  =  0 

3t  M  i,j  J  y  ,1  Re  \  i,j j  3  j,ji ) 


II  _  Hi  3R  + 

dt  Y  3t 


PT  .u. 
,1  i 


Y-l  ,  y(y-1)  /  2 

— 1 —  p  . u .+  -  —  1  —  u.  .  u.  .  -  u. 

Y  ,ii  Re  \3  1,1  j ,j  i,j 


-  u.  . 
J.i 


p  =  pT 


(3) 

(A) 


where  the  commas  denote  partial  derivatives,  the  repeated  indices  imply 


summing  with  the  range  of  index  i  or  j  being  3.  We  use  index  notation 
in  preference  to  vector  symbols  in  order  to  facilitate  clarification  in¬ 
volved  in  integration  by  parts  and  computer  coding  in  our  later  dis¬ 
cussions.  The  nondiraensional  quantities  are  defined  as • 
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The  presence  of  viscous  terms  is  intended  for  development  of  shear 
layers  close  to  the  burning  surface  although  the  effect  of  viscosity 
is  negligible  in  the  flow  domain  with  high  Reynolds  numbers. 
Substituting  Eqs.  (1)  and  (4)  into  Eq.  (3),  we  obtain 


3t 


v  2  /  2  \ 

+  Ypu  •  -+P.U.+—  -ru.  .u.  ,-u.  .u.  ,-u.  .u.  .  \ 
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=  0  (5) 

To  obtain  the  acoustic  equation  we  take  the  spatial  derivative  of 
Eq.  (2)  and  subtract  the  results  from  the  time  derivative  of  Eq .  (5). 


32d  £ 

3  "P-  p 


p.u  ' '  h  (»,!  "i)  uj),i] 
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(6) 


At  this  point  we  introduce  the  perturbation  expansion  to  0  (£ )  for 


the  velocity  u^,  pressure  p,  and  density  p  in  the  form. 


u.  =u.  +  £  fu.  +  u  ^ 
l  i  V  i  1/ 


p  =  1  +  ep' 

p  =  1  +  ep ' 


(7) 

(8) 
(9) 


where  the  bar,  asterisk,  and  prime  denote  the  mean  flow,  vortical  fluc¬ 
tuation,  and  acoustic  fluctuation,  respectively.  Expanding  Eq.  (6)  in 
terms  of  the  perturbation  variables  and  collecting  the  0  (c)  terms,  we 
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arrive  at  the  expression: 
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Considering  the  vorticity,  defined  as, 

£ijk  \,j  ^i  ’  £ijkUk,j-  ^i  ’  ^i  ~  ^i  +  ^i 
the  acoustic  equation,  Eq.  (10),  is  recast  in  the  form 

P%,  -  IV  ,  u 


~  3t2  h 


where 
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with  being  the  permutation  symbol.  The  vorticity  transport  equation 

1 J  k 

is  obtained  by  taking  a  curl  of  the  vector  form  of  Eq.  (2),  and  collect¬ 
ing  the  terms  of  0  (e), 

jV 

J  r_  *  *  ,\_T  1* 

-r- —  e. ..  e,  u  £  +  (u+u)£  .-  —  5.  ..*0  (14) 

3t  ljk  kmn  [_  m  n  V  m  m  /  n_|,j  Re  i,jj 

Equations  (11),  (12),  and  (14)  represent  the  most  general 


vorticity-coupled  acoustic  problem  in  which  the  viscous  action  of  the 
fluid  is  fully  taken  into  account.  Here  the  main  objective  is  to 
determine  the  growth  rate  given  by  the  stability  integral.  This  sub¬ 
ject  is  discussed  in  the  following  section. 

3.  THE  STABILITY  INTEGRAL 

3.1  Vorticity-Coupled  Acoustic  Instability 

To  begin  we  must  first  consider  the  boundary  conditions  at  the 
burning  surface.  In  the  direction  normal  to  the  surface,  the  momentum 
equation,  Eq.  (2),  takes  the  form,  in  terms  of  0  (e) , 

-"'i  ”i  '  £  05) 

with 


f  ■  ui  +( 

u.  u .  \  .  +  (u.  u:)  .-£...  (u.  E,*  +  u*  L  + 
J  JAi  V  3  3/,i  xjk  \  j  sk  j 

uj  0 

i 

J 

-  h  If*  *  ui ) 

..  +  ^  (u*  +  ur)  ..jin. 

OJ  3  \  3  J  J  |  i 

(16) 

i 

The  oscillatory 

motion  of  the  acoustic  media  is  modeled  by 

i 

--  „  ikt 

p  =  pe 

- 

(17a) 

-  /v-  ikt 

u.=  u.e  , 

1  X 

*  ikt 

u .  =  u.e 
i  i 

(17b) 

r*  r*  ikt 

V  cie 

(17c) 

, 

where  k  is  the  complex  dimensionless  frequency  given  by 

k  =  oj  -  ia  (18) 

Here,  the  imaginary  part  is  known  as  the  growth  rate.  The  instability 

cc  t 

of  acoustic  pressure  is  signified  by  the  growth  of  a  proportional  to  e 

7  ' 
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Substituting  Eq.  (17)  into  Eqs.(12)and  (15)  yields,  respectively, 

h  -  iku  £  +  ikYu  p  -  y|(u  u*)  +  (u.  uT)  .  -e.  (u.  t*  +  u*  L 

1  i.1  [_v  J  ]  '  j  J/,ii  ljk  V  j  k  3  k 

+  Uj  + fe[(G* +  +  k  (s* +  g;  ),ju] 

-  u.  .  ( u.  +  uT  )  .  (19) 

J  j1  '  3  J  '  ,i  K 


f  =  Y  )  ikuC  +  (  u.  u.  )  +  (  u.  uT  )  .  u.  £*  +  u.  £. 

{  i  v  3  J-'.i  '  J  i]kV  j  '’k  j  sk 

+  5-  U-±T«  +  «) +4(s*  +  g:)  ,.]L.  c 

j  ^k'  Re(_\  i  i',jj  3  \  j  f  i  ^ 

The  foregoing  process  leads  to  a  nonhoraogeneous  Helmholtz  equation 

P,ii  +  k2p  =  h  «  (I 

subject  to  the  boundary  condition 

pVi  - 4  (2 

For  the  solution  of  Eq.  (21)  we  make  use  of  the  Green's  function 
[8] .  It  can  easily  be  shown  that 

(k’-k^Ej.  /fipNdB+  /  f  P„  dr  (2 


where  the  unperturbed  mode  shape  p^  and  the  wave  number  k^  are  determined 
from  the  classical  acoustic  problem, 

fs.il  +  k’f»  •  0  <24) 

f»,i\  •  0  <25) 

Note  that  is  given  by  the  integral 

E>  .  /flf;  dS!  (26) 
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It  should  be  noted  that  a  correct  integration  of  the  *.  ght-hand  side 
of  Eq.  (23)  is  the  most  crucial  aspect  of  the  present  study.  That  is, 
integration  by  parts  must  be  carried  out  until  all  Neumann  boundary  con- 

A 

ditions  are  brought  to  the  surface.  Since  h  in  Eq .  (19)  contains  the 
terms  from  the  momentum  equation  differentiated  once, it  is  now  necessary 
that  they  be  integrated  by  parts  twice  in  order  to  produce  the  Neumann 


boundary  conditions.  For  example,  consider  a  single  term  taken  from 
Eq.  (19)  substitued  in  Eq .  (23)  and  integrated  by  parts: 


(  (-(ITu.)..  p  dft  =f(u'u.).  n.  p  dT  -f(u'u.).  p  . dft 
]  J/,ii  N  JT\  2  1  N  2  J,i  N,l 


(27a) 


Integrating  by  parts  now  the  second  term  of  Eq.  (27a)yielas 


■  /  ( u'u.  )  .  p  .  dft  =  -f  u  .n.  p  .dT  +  fuT  u.  p  .  .  dfi 
jr;  J  jU  N,1  Jr  J  2  1  FN,1  Jn  2  1  N,n 


T27b) 


In  order  to  demonstrate  the  advantage  of  tensor  notations  over  vector 
notations  in  the  process  of  integrations  by  parts,  we  consider  a  following 
example  (eighth  term  of  Eq .  (19)  substituted  in  Eq .  (23). 


A  A  -  A 

I'  u.  .  .  ,P  dfi  =  f  u.  .  ,n.P„  dT  -  u.  ,PM  .  dS 
i,jji  N  J  i,jj  l  N  Ja  i,j:  N , l 

l  T  u 

=  f  u  ,,n,P  dT  -  f  u.  .n.P  dF  +  f  u.  P 
jr  i  N  Jr  i,j  2  N,x  Ja  1,3  N,i; 


(28a) 


If  we  pursue  the  above  process  in  vector  notations,  we  will  encounter 
difficulty  in  obtaining  the  last  term  of  Eq.  (28a).  To  see  this  we  begin 


with 


f  V  •  V2u  P  dQ  = 

h  ~  ~  N 


’  v2  «  pNdr 


f  r2u  •  V  P„  dQ 

h  ~  ~  N 


=  j  n  •  V2  u  P^  dF  -  J  (n  •  V)  u 


!  PN  dF  + 


f  (  ?  )  dQ 

J  o 

‘  (28b) 
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In  integrating  by  parts  /  ^2u  •  V  P  dfi  ,  only  the  boundary  term  can 

'n  '  -  N 

be  obtained  in  vector  notations.  To  obtain  the  domain  term,  the  quantity 

A 

V2u  .*  ^T>n  must  be  expanded  into  individual  terms.  Other  than  in  special 
cases  such  as  this,  however,  the  vector  notations  can  be  used  to  perform 
integration  by  parts.  But  care  must  be  exercised  to  ensure  correct 
product  operations. 


An  evaluation  of  the  rest  of  the  integrals  in  Eq .  (23)  leads  to 
lk  "kN)EN=:CYkN  4  UiPNnidr  +  lkN(Y+1)  UiPNnidr  “  lkN  4  “i.^N^ 

-ikN(2Y+l)  VNPN,idn  +  yl  “j(Gj+Gj)PN,inidr  ‘  yi  "j(Gj+Qj)PN,ii‘ 


- 

,  £  .  . .  (  u .  £,  +  u .  L  +  u 

1  ijk  \  j  k  jTc 

:)  PN,idn 

-  JL 

Re 

[  ff''*  4.  ■'') 
r  L'ui  +  V.j  pN,inj  + 

(u.  +  uT)  .pN  .n. 
3  '  j  j /, j  N, l  l 

4  [(GI  +  Gi),jpN,ij  +i 

(u.  +  u')  ,p  .  Idft 

2iky(Y-l)f  f2  -  /.*  \ 

Re  Jn  i-3  Ui,i(Uj  +  Uj),j 

-  u.  .  (u*  +  uA  . 

\  J  -  J/,1 

u . 

-  3  * 

i(“j  +  aJr),l]?»dfi 

A  A 

Note  that,  in  view  of  our  choice  made  in  Eq.  (25),  we  set  p  **  PN  and 

k  =  in  Eqs .  (19)  and  (20).  Squaring  both  sides  of  Eq .  (18)  and  in 

view  of  o  =  a)  =  k  at  a  =  0  for  neutral  stability,  we  obtain 
N  N 

k2-  k^  ”  -  i2otk^  +  a2  (30) 
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Comparing  Eq.  (30)  with  Eq.  (29),  we  may  solve  for  the  growth  rate  a 
(a^<<  1 2al^|  )  by  equating  the  imaginary  parts. 


where  the  superscripts  (R)  and  (1)  refer  to  the  real  and  imaginary  parts, 
respectively,  and  the  various  terms  are  defined  as: 

(A)  Surface  combustion 

(B)  Surface  convection 

(C)  Surface  viscous  damping 

(D)  Combustion  into  domain 

(E)  Convection  into  domain 

(F)  Momentum  viscous  damping 

(G)  Dissipative  energy 


1L 


It  should  be  noted  that  the  mean  flow  vorticity  £^ ,  the  vortical 
fluctuation  £^,  and  the  vortical  component  of  the  velocity  u^  must  be 
known  in  order  to  evaluate  the  integral,  which  we  shall  discuss  in 
Section  3.2.  It  is  possible,  however,  to  evaluate  the  stability  integral 
in  an  alternative  form  without  explicit  values  of  the  vorticity.  To  this 
end,  we  proceed  with  the  convective  term  (u  *  V)u  instead  of  V (y  u2)-u  x  £. 
Thus  we  return  to  Eq.  (23)  and  examine  the  integral  of  the  form 

V  •  (u  •  V)  u  pNdfi  =  (Ui,j  Uj),i  ^Ndn  (32) 

Integrating  Eq.  (32)  by  parts  twice,  we  obtain  the  growth  rate  similar 
to  Eq.  (31)  except  that  the  terms  designated  as  (B)  and  (E)  are  replaced 
by 

°<B)'-2iT  |“i(aj+a  j)  +  +  al)a>  "jlfN,injdr  (33a) 

1  f  y  )f-  M  *  aA(i)  M*  ,  ~\(D- 

a  =  -  -r— y  1^  -  ; —  <  u.lu.  +  u.  (  .  +\u.  +  u.  I  u,  .  p.T  . 

(E)  2E^  JSl  iV  j  J', J  '  1  1 '  j,j  i  N,i 

+  l  Gi  (°j  +  Gj)U)  +  (°i  +  Gi)  "j  1  ^N,ij[d^  (33b) 


A  careful  examination  of  the  integrals  above  is  in  order.  Recal  that  the 

counterparts  of  a,,,.  and  a  in  Eq .  (31)  were  obtained  from 
<B)  (E) 

|  V  .  [7(i»’)-  »  x  C)  P>-/  VO  (34) 

n  n 

The  vorticity  can  be  easily  removed  from  Eq .  (34)  by  setting  £  =  0.  This  can 

~  * 

not  be  done  in  the  case  of  Eq .  (32).  The  vortical  component  u  can  be  set  equal 
to  zero  in  both  Eqs .  (32)  and  (34)  with  £  =  0  in  Eq .  (34),  but  this  will  not 
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make  both  equations  equivalent.  In  other  words,  the  effect  of  irrota- 
tionality  cannot  be  assured  from  Eq.  (32)  simply  by  setting  u*  *  0,  and 
we  are  aware  that 

[(u  +  eu")  •  V] (u  +  ejT)  +  (u  +  eG")  •  (u  +  eG')] 

If  the  viscous  effect  is  neglected,  then  the  terms  designated  by 
(C) ,  (F) ,  and  (G)  will  be  eliminated.  Thus 

“  -  -  2Ej  [(A)  +  (B)  +  (D)  +  (E)]  (35; 


Evaluation  of  the  integrals  given  by  Eqs.  (31)  and  (35)  requires' 
adequate  analytical  or  numerical  models  for  the  mean  flow  and  vortical 
and  acoustic  components  of  the  velocity.  These  topics  will  be  discussed 
in  Section  5. 

If  we  assume  =  0,  u*  =  0,  and  Re  ->•  »  in  Eq.  (31),  we  obtain 

(ct  )  =  -  —  r  f  r  ~  * 

a  Ci=0  2En  \y  PN  n.  +  (Y  +  1)  u.P^n. 


Y  "*  (IJ  “1  r  r  -  *  ~  /s 

+  s,  ”i  “j  .V  "1  ]dr  +  /n['ui,i  ?i  -  <2Y  -  »  :1  V>,1 


'(I) 


rf-  u.u;(I)  p 


».u3d!!] 


(36) 


Here  the  velocity  normal  to  the  boundary  surface  is  expressed  in  terms 
of  the  admittance,  A=A^  +  iA^^  and  the  mean  flow  Mach  number  M 
such  that 


uTn.^fuT^+i  uT^\n.  =  M  A. ,P„/  y  (37) 

1  1  V  i  1/1  N 

whereas  the  acoustic  fluctuation  velocity  in  the  domain  is  given  by 


U. 

X 


V  pNpi 


(38) 
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Invoking  the  relationships  of  Eqs. 
(incompressible  flow) ,  one  obtains 


(24)  and  (25)  and 


P..  +  u.  P. 


2  n  •  1 
N  ij 


df 


setting  u. 

i  1 


=  0 


(39) 


which  is  the  familiar  expression  as  derived  by  Culick  [see  Eq.  (2.29)  of 
Ref.  4  neglecting  particle  distributions].  This  is  the  simplest  form 
representing  the  possible  unstable  motions  in  a  solid  propellant  rocket 
motor.  It  has  the  ingredients  such  as  velocity  coupling  (burning  rate 
changes  in  response  to  velocity  fluctuations  parallel  to  the  surface)  and 
pressure  coupling  (burning  rate  changes  in  response  to  pressure  fluctua¬ 
tions  on  the  surface) . 

Culick,  [3,4]  further  defines  the  "flow  turning"  as  an  averaged  appro¬ 
ximation  to  viscous  effects  which  appeared  in  one-dimensional  analysis  but 
did  not  arise  in  the  three-dimensional  inviscid  flow.  This  claim  has  also 
been  substantiated  by  the  present  authors  as  shown  in  Eq .  (39). 

Let  us  now  investigate  the  terms  which  result  from  Eq .  (33a, b).  The 
stability  integral  for  the  case  of  Re  -*■  00  assumes  the  form 


a 

a 


-  2Ej  I  4  lY  fii<R)Vi  +  (Y+1)  "i  pn  rx 

L  y  (  -  x-(I)  .  -/(I)-  V  n  I 
+  ,  l  U .  u .  +  u .  u .  )  p.T  .  n .  dr 

k„  \  i  J  i  1  J 


u.  . 

1,1 


(2yt-l) 


UiPNPN,i 


f(;. a:(»+ .) 

b* '  1  ij  1 


N,i 


a-(I) 

J 


(40) 
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Now,  using  the  relation  given  by  Eq.  (38)  ,  we  have 

2^1  Jr  K'(I\  "i  +  ", 

+  4  ("i  VjVi  +  "j  YiVl’  "j  ] dr 

+  I  *  ui,iPN  ~  (2”V+I)  ui^N^N , i  “  (  ui  jj  ^N,i 

+  "j,j  ViVl  +  “i  +  “3  ^»,A,ij)l 

For  u.  .  =  0,  and  using  Eqs.  (24)  and  (25),  we  obtain 

1,1 

a  *  -  -r^r  F  f  (yu'^ n,PM  +  u  P*  n  +  A-  u.P  .?  ,  n.  ) 
a  2E^  j'j,  I'  i  i  N  i  N  i  k‘  j  N,i  N,i  J  ) 


(A-2) 


^N  fn  (Ul^N*^N»id  +  UjPN,iPN,i1)JC 


According  to  Culick  [3,4],  the  first  two  boundary  terms  are  defined  as 
combustion  (at  surface)  and  the  last  boundary  term  is  known  as  the  mean 
flow/acoustic  interaction  or  "flow  turning".  "Note  that  the  "flow  turning" 
term  has  appeared  from  the  integration  of  Eq .  (32)  by  parts  "twice".  The 
boundary  term  which  results  from  the  first  integration  cancels  with  a 
term  from  a  boundary  integral  of  Eq.  (23).  Upon  integrating  by  parts 
"once  more",  we  obtained  the  boundary  terms  in  Eq.  (33a)  which  finally  led 
to  the  "flow  turning"  term.  Recall  that  it  arises  from  the  convective 
term  of  the  momentum  equation.  The  same  term,  designated  as  "B"  in  Eq .  (31) 
was  called  "surface  convection".  Conversely,  the  domain  terms  arising  from 


*  .  w* .JSg. 


the  convective  terms  of  the  momentum  equation,  designated  as  "D",  were 
called  "Convection  into  domain".  These  convection  terms  correspond  to  the 
popular  definition  of  the  "flow  turning"  as  a  consequence  of  injecting 
mass  into  the  acoustic  wave,  which  must  supply  acoustic  energy  to  the 
mass  in  order  that  it  be  brought  up  to  the  local  velocity  [3,4].  Such 
action  must  combine  both  surface  and  domain  integrals,  linking  the  surface 
activities  into  the  domain. 


3.2  Vorticity  Stability  Integral 

Recall  that  the  coupling  of  vorticity  with  an  acoustic  field  given  by 
Eq.  (31)  is  yet  to  be  evaluated  from  Che  acceptable  velocity  profile  in 
the  shear  layer.  The  influence  of  direct  acoustic  feedback  and  presence 
of  unstable  shear  flow  are  to  be  taken  into  account. 

One  option  for  the  solution  is  to  establish  an  independent  stability 
integral  for  vorticity.  To  this  end,  we  return  to  Eqs.  (14),  (17b),  and 


(17c),  and  write 

j[  yN*  ^  A  ys 

Re  Ci.jj  +  =  Hi 


where 


Hi  eijkEkmn  [uj*n  +  (0*  +  \  I  f j 
subject  to  the  boundary  condition 


_i 

Rti' 


(43) 


(44) 


L6 


It  is  evident  that  Eq.  (43/  represents  diffusive  waves  due  to  production 
of  vortex  streets,  and  that  the  state  of  instability  is  once  again 
determined  from  the  Green's  function  technique  [8] .  It  can  easily  be 
shown  that 


( k  -  s)  -  4  kKda + 4  sfdr 


where 


=  L  ^i  ^i  ^  =  h  £<4V  £  i  }.m  4G!NJ.  dQ 


Jn  ijk  i&m  k,j  m,£ 


The  unperturbed  mode  shapes  and  the  wave  number  k^  are  calculated 
from  the  homogeneous  equations 

t  ?i.”  j  +  k  ?rB-  °  (< 


The  vorticity  growth  constant  a can  be  derived  by  the  relation  (18) 


and  (-4  5), 


f  £  e  +  g*N|  )£*> 

ijk  kmn\  m  n  m  n/  x,j 


Note  that  the  term  associated  with  acoustic  fluctuation  velocity  u^  in 
Eq.  (4&)  drops  out  as  a  result  of  squaring  the  imaginary  number  which 
appears  also  in  Eq.  (38),  thus  making  this  term  a  real  part.  If  desired 

for  a  computational  standpoint,  we  may  replace  the  vorticity  in  terms  of 


velocity. 


p*N  „  0*N 

^i  ijk\,j 


(49a) 


Ci  -  'ijkVj 
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(49b) 


Evaluation  of  Eq.  (46)  or  Eq .  (48)  is  now  straightforward  once  the 
mean  flow  u^  and  fluctuating  velocity  u^  or  the  vorticity  can  be 

numerically  determined  using  the  finite  element  method.  The  hyperbolic 
tangent  velocity  profile  for  the  mean  flow  and  shear  layers  [9-12]  may 

be  used  for  one-dimensional  integration.  I 

The  growth  rate  for  the  vortex  shedding  cc  ^  represents  an  independent  , 

vorticity  instability  equivalent  to  the  solution  of  Orr-Summerfeld  equation. 

i 

When  vortex  shedding  frequency  approximates  the  classical  acoustic  frequency, 
periodic  flow  separations  can  produce  significant  pressure  oscillations. 

! 

To  this  end  we  perform  the  eigenvalue  analyses  for  both  Eqs.  (24)  and  (46)  j 

and  we  enter  the  results  into  Eq .  (31). 


4.  FINITE  ELEMENT  ANALYSIS  ; 

I 

I 

j 

The  use  of  finite  elements  in  fluid  mechanics  problems  has  in-  j 

creased  significantly  in  recent  years.  An  application  of  this  method 
to  combustion  instability  was  first  studied  in  [7],  Our  current  effort 

is  directed  toward  extending  the  work  of  [7]  to  incorporate  the  new  terms  , 

in  the  stability  integral  as  described  in  Section  3.  j 

To  begin  we  return  to  a  classical  acoustic  problem  characterized  by  j 

Eq.  (24)  and  Eq.  (25).  The  Galerkin  finite  element  equation  takes  the 

j 

i 

f 
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form  [13] , 


4(p,ii  +  k2  pk<^  =  ° 


(50) 


where  is  the  test  function  which  is  set  equal  to  the  trial  (basis) 


function  such  that 

P  (  ?,*)  =  (5)  p 


(51) 


Here  x  and  t  denote  spatial  and  temporal  variations,  respectively,  the 
subscript  a  representing  the  global  node  number  with  a  =  1,2,  -  -  -  ,  n, 
n  being  the  total  number  of  global  nodes. 

It  follows  from  Eq .  (50)  that  the  finite  element  eigenvalue 
equation  is  of  the  form 


aB 


-  k  B 


aB 


where 


aB 


4  *a,i  ^8,i 


df2 


aB 


l 


o  *  *0  dft 

ft  a  0 


(52) 


(53) 


The  normal  modes  pN  required  for  the  stability  integral  can  then  be 
determined  from  Eq.  (51).  Non-axisyrametric  geometries  for  slots, 
segments  or  irregular  flow  domain  are  modeled  routinely  using  the  Fourier 
series  expansion  [7] . 

Similarly,  the  vorticity  transport  Eq.  (  4G)  is  cast  in  the  form 

1=  0  (54) 


Aik  -  k  Bik 
aB  aB 


where 


Aik  =  A.  f  $  $  5  dft  ,  Bik  =  4  Vb' 

aB  Be  h  a,j  B,j  ik  * 


6.,  dft 

ik 


(55) 


The  well-known  QR  algorithm  [14]  is  invoked  for  the  solution  of  eigenvalues 

*  N  , 

and  eigenvectors.  As  a  result,  the  normal  modes  or  are  calculated 
and  substituted  into  the  stability  integrals,  Eq.  (31)  and  Eqr(48). 
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Once  the  normal  modes  are  determined  either  numerically  or 
anlaytically ,  the  evaluation  of  stability  integrals  by  finite  elements 
is  most  appealing.  Highly  nonlinear,  complicated  functions  can  be 
integrated  quite  accurately  via  Gaussian  quadrature, using  the  iso¬ 
parametric  finite  elements  [13]. 

For  example,  the  stability  integral  is  written  in  the  form 


G  <£,  n,  £  )  dn  dC 


(56) 


This  leads  to  the  Gaussian  quadrature  process. 


m 

E 

m 

E 

m 

w.w  wk  g«  n  ck) 

(57) 

1=1 

j=i 

k=l 

where  the  weighting  functions  w.  w.  w  and  abscissae  £.  n.  C  are 

r »  J , 

chosen  for  an  adequate  number  of  Gaussian  points . 


5.  APPLICATIONS 


5.1  General 

The  stability  integrals  given  in  Section  3  and  the  finite  element  ' 
equations  as  shown  in  Eq.  (57)  will  be  discussed.  We  begin  with  boundary 
conditions  defined  on  the  burning  surface  (1^)  and  the  nozzle  entrance 
(r„)  as  represented  by  the  admittance  functions: 


Ui  ni 


^  *  "  %  p7y 


on  r. 


(58) 


_Ui  ni 

Mn  P7  Y 


on  r 


(59) 
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where  and  Mn  are  the  mean  flow  Mach  numbers  on  P^  and  T  ,  respectively. 

At  this  point,  it  is  important  "o  realize  that  the  different  forms 
of  convective  terms  in  the  momentum  equation  would  give  different  results 
for  combustion  instability  phenomena. 

Case  1:  u2)  “  u  x  £  (60) 

Case  2:  ^u  •  vju  (61) 

The  difference  is  particularly  significant  when  dominated  by  vortex 
sheddings.  To  see  this,  we  rewrite  Eq.  (23)  in  terms  of  admittance 
functions  for  the  Cases  1  and  2,  respectively, 

Case  1: 


a'%  +  aH 


(62) 


where 


a  =  - 


3 


a(r)m 
n  n 

ar„  *  4  fji*« 

n 

1  -  X 

PN,i  PN 

n.  +  4  PM  .P.T  .  .n.)dr 

J  3  N,i  N,jj  i J 

+  ir  -  (2y+i)  si  %  pN>i  -  4  l^j  "N,j  pN,ii 

+  £ijk  ^k  PN,i  PN,j  J  +  k*Re  (  PN,ijPN,ij  +  3  PN,iiPN,jj  ) 


(63a) 


“H  =  '  2E3 


f  j  u  u*^P  n - —  /  u*  p  n  +  —  U*(I>P  n  'jidf 

/Sj  J  2  N,i  1  kNRe  ■  ^  N,i  2  3  N,i 


+  f  I'  ^  [u.u*(I)p^.  +e..  /u.^Vu*^  ?  U  I 

Jnl  ^  j  J  N,i!  i]k\  j^k  3  sk/rN,iJ 


J  a 


.  X/G*(I)  fi  +  A  n*(I)  -  Y^-I^T2  "  ~*(R) 

kVRe  1,3  Pn»13  3  30  PN,ii/  Re  \3  Ui,i  Uj , j 

-  r*(R)  -  **(R)  \  „  I  ,0 

i, j  J,i  J , i  3,1  /  KN  I 


(63b) 
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> 


Case  2 : 


°*  '  '  2*i 


-  ^  «b  L  dIb  +  AbW  fin  lr  K  drn 


(A-l) 


-  (T  +  l)  Jr  ^  "i  dr  +  4  P»,i  V  “lFM,iP».jn3 


dr 


(A-2) 


(B-l) 


(B-2) 


*  1 

&  j  (F»,i  S»,ij  "j  +  J  P»,jjVi"i)dr  +  f  I'  =i,i  ^  -  (2Y+1)  :1  %  Vi 

r 


(C)  (D-l) 

~  kj  ^  Ui  ^N,jj  **N,i  +  Uj,j  **N,i  ^N,i  +  uj  ^N,j  ^N,ij 


(D-2) 


(E-l) 


(E-2) 


(E-3) 


+  uj  ^N,i  ^N,ij)  +  k^Re  (  ^N,ij  %,ij  +  3  ^N,ii  ^N.jj  ) 


dQ 


(E-4) 


(F-l) 


(64a) 


“«  ‘  ~  2EJ  [^{i^^uiuj  +  u  '  ’“PVi"  kBRetUi,jP»,inj 
+  dr  +  jj  '  ^  (“i  V?  %,i  +  =j,j  ai”<r,pN>1 

+  "i  S*a>5».ij  +  =3  “I'1’  *<«.«)  +Ve(  ai!j)p»,i3 

-*(R)  \  a 
-  u.  .  u.  .  pM 
J,i  J  > 1  /  N 


(64b) 


A  « 

Here  cl,  denotes  the  coupling  of  vortical  velocity  fluctuations  u^  and  the 


“h 

local  acoustic  pressure  p 


N 
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5.2  Vorticity-Uncoupled  Acoustic  Instability 

We  consider  the  axisymetric  conical  geometry  (Figure  la)  with  a 
transitional  angle  (p  and  a  port  diameter  D  which  diverges  into  2D  at  the 
nozzle  entrance,  the  transition  beginning  at z  =L/3.  It  is  assumed  that 
aM^  =  10  in/sec,  y=  1.2,  and  V  =  0.1  in2/sec.  It  is  further  assumed 
that  the  mean  flow  velocity  is  given  by  [16] , 


(65) 


where  R  is  the  radius  as  a  function  of  z  and  e  and  e  are  unit  vectors 

r  z 

in  the  radial  and  axial  directions,  respectively. 

The  classical  acoustic  modes  pN  are  given  by  the  formulas 

PN  -  cos  (kzz)  cos  (m0)  Jm  (k^r)  (66) 

with 

^  *  k£  +  (67) 

where  m  *  0,  1,  2, - »  k£  =  with  *  =  °»  1»  2> - ,  and 

k  are  the  roots  of 
mn 

Jm  (k»r)  =  0  at  r  =  R  (z) 

To  evaluate  in  Eq.  (64a)j the  vorticity-uncoupled  growth  rate,  we' 
make  use  of  Eqs,  (65-67),  rewrite  Eq.  (6Aa)  in  terms  of  finite  elements,  ' 
and  transform  the  integrand  into  the  form  of  Eq.  (57).  To  demonstrate 
the  accuracy  of  calculations  we  choose  the  results  of  terms  given  by 
(A-2)  and  (B— 1)  which  are  shown  in  Table  1.  Here  we  used  4x4  Gaussian 
points,  <p  =  0,  L/D  =  1  ,  A^  *  1,  A^  M^/M^  =  A  comparison  with 

analytical  calculations  indicates  an  excellent  agreement. 


24 


Mode 

a(A-2) 

a(B-l) 

(£mn) 

(ESHSsSI 

Present  Study 

Present  Study 

(100) 

-4.400 

-4.400 

2.000 

2.000 

(200) 

-4.400 

-4.400 

2.000 

2.000 

(010) 

-6.253 

-6.253 

0.849 

0.848 

(110) 

-6.253 

-6.253 

1.683 

1.683 

Table  1.  Acoustic  Instability  Growth  Rate  -  Integrals  (A-2)  and 

(B-l)  in  Eq.  (64a). 

Figure  2  shows  the  acoustic  instability  growth  rate  ag  versus  the 
transition  angle  4>  for  several  inodes.  It  is  interesting  to  note  that,  . 
for  both  Case  1  and  Case  2,  the  tangential  mode  tends  to  be  more  unstable  in 
comparison  with  other  modes,  and  that  instability  increases  with  larger 
transition  angles.  However,  for  a  given  mode.  Case  1  exhibits  more  in¬ 
stability  than  Case  2.  To  the  best  of  our  knowledge  this  has  not  been 
brought  to  the  attention  of  the  combustion  community  [17].  Furthermore, 
it  has  been  demonstrated  in  Table  2  that  the  effect  of  viscosity  (C+F) 
is  negligible  in  the  cylinder  of  uniform  diameter  (no  flow  separation). 

It  should  be  noted  that  an  alternative  approach  is  to  perform  tTie 
eigenvalue  analysis  for  Eq.  (52).  The  acoustic  modes  pN  are  then  cal¬ 
culated  and  used  in  the  evaluation  of  stability  integral.  For  irregular 
geometries  (non-axisymmetric)  such  as  occur  in  the  star-shaped  pro¬ 
pellants  we  must  resort  to  the  eigenvalue  analysis.  Numerical  results 

using  the* procedure  as  employed  in  [7]  are  forthcoming  in  a  subsequent 


paper . 


^\Integrals 

Mode's. 

(£mn)  n. 

A 

3 

C 

D 

E 

F 

(001) 

Case  1 

-.0762 

0.0 

-.0000 

.6793 

-  .1998 

-.0421 

Case  2 

-.0762 

-3.7400 

-.0000 

.6793 

-  .5994 

(010) 

R 

5.3288 

0.0 

0.0024 

-2.3040 

0.6777 

-.1102 

0.0 

0.0024 

glllllll 

BSBI 

RB 

(100) 

Case  1 

-7.0292 

0.0 

0.0061 

6.4790 

-1.9060 

-.0221 

Case  2 

-7.0292 

2.9130 

6.4790 

-7.561 

-.0221 

(110) 

Case  1 

-1.4362 

0.0 

0.0035 

2.9720 

-  .3740 

-.1144 

Case  2 

mm 

Igjggfcg™ 

(200) 

Case  1 

-4.8580 

0.0 

0.0060 

4.1420 

-1.2180 

-.0225 

Case  2 

BojXuB 

BBSil 

- .0225 

Table  2 

.  Growth 

*  L 

Rates  for  —  =  4 

,  <j>  =  ~ 

*  ”  4  > 

Eq.  64a 

5.3  Vorticity-Coupled  Instability 

A  complete  three-dimensional  analysis  for  the  vorticity-coupled  in¬ 
stability  must  follow  the  eigenvalue  analysis  as  required  by  Eq.  (54). 
The  vortical  modes  and  disturbances  thus  calculated  may  be  sub¬ 
stituted  into  oty  in  Eq.  (63b)  or  Eq.  (64b).  Similar  calculations  for 
Eq.  (47)  can  be  carried  out  for  the  case  of  an  acoustic-uncoupled 
vorticity  instability. 

In  this  paper,  however,  it  is  our  plan  not  to  be  involved  in  the 
eigenvalue  analysis  but  to  show  an  analytical  approach  using  the  hyper¬ 
bolic  tangent  velocity  profile  for  a  shear  layer  [10-12]  .  Here  we 
assume  the  stream  function  ip*  (r,z,t)  representing  a  single  oscillation 
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Transition  Angle  <£ 

Figure  2  The  Acoustic  Growth  Rates  Vs._Transition  Angles  for 
Different  Modes,  L/D  =  A  and  aM^  *  10  in/sec. 
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of  a  disturbance  to  be  of  the  form  [9] . 


A  a  k  ilct 

ip  (r,z,t)  =  Ip  (r,  z)  e 


where  k  is  the  complex  dimensionless  frequency  defined  in  Eq.  (18).  We 
further  assume  that  the  disturbance  amplitude  of  the  stream  function 

A* 

ip  (r,z)  has  the  same  form  as  in  [10] 

£  (r,z)  =  e  e1®2  (69) 

where  y  and  8  are  real  quantities  representing  the  various  mode  shapes 
which  may  depend  on  the  geometry,  and  y"  and  z'  are  the  coordinates 
normalized  with  respect  to  the  boundary  layer  thickness  0q  at  the 
separation  point  ( z "  =  0)  and  measured  from  the  inflection  point  as 
shown  in  Figure  la.  That  is,  the  disturbance  effect  of  vortex  shedding 
on  the  acoustic  field  is  assumed  to  begin  at  the  critical  point  (z'~  0) 


and  vanish  downstream  before  entering  the  nozzle.  Figure  lb  shows  the 


disturbance  function  tp  ^  )  coupled  with  the  acoustic  field  pN  for  a 

/v*  ** 

longitudinal  mode.  The  velocity  disturbances  uz  and  uf  for  an  incompress¬ 


ible  axisymmetric  flow  are  given  by 

*k 

l  dip 

U  =  —  •*"  ••  # 

z  r  3  r 


i  dip 
ur  "  _  r  3z 


The  mean  flow  representing  a  hyperbolic  tangent  velocity  profile  [10-12] 
assumes  the  form 


u  ■  4  (1  +  tanh  y') 
z  2 


The  shear  layer  thickness  0  is  found  experimentally  to  be  a  function  of 

o  • 
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(72) 


Reynolds  numbers  and  the  jet  diameters, 

e  --£5l_ 


✓_  (S) 
Re 


where 

Re(S)  =  U  D/v  (73) 

o 

with  Uq  being  the  free  stream  velocity,  and  cq-0.62  for  cylindrical 
geometries  [12] . 

For  various  vortical  mode  shape  combinations  of  y  and  B  the  growth 
rates  and  are  calculated  for  the  two  first  longitudinal  modes 
(Table  3).  Here  the  same  physical  constants  are  used  as  in  Figure  1 
with  <{>=tt/2.  Note  that  is  independent  of  the  acoustic  mode  and  0^=0 
for  any  tangential  modes.  The  results  indicate  that  is  strongly 
dependent  on  6  while  slightly  dependent  on  y  for  the  two  longitudinal 
modes,  and  that  the  higher  mode  is  more  stable  than  the  lower  one.  Note 
also  that  the  trend  of  is  similar  to  that  of  but  has  the  largest 
instability  occurring  at  6  »  0.6,  and  neutral  stability  at-3^1  for  all  y's. 

For  the  first  longitudinal  mode,  the  values  of  aH  are  in  the  range  of 
the  case  investigated  in  [9] . 


6.  CONCLUSIONS 

A  significant  improvement  over  the  current  practice  of  combustion 
stability  calculations  is  achieved.  A  three-dimensional  formulation  of 
stability  integral  for  the  coupled  acoustic-vortex  fluctuations  leads  to 
new  integral  terms  accounting  for  the  various  phenomena  previously 
neglected . 
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One  of  the  terms  in  three-idimensional  combustion  instability 
integral,  identified  as  "flow-turning”,  was  shown  to  be  the  result  of 
integrating  by  parts,  twice,  one  of  the  convective  terms  of  the 
momentum  equation.  The  most  crucial  factor  is  the  correct  procedure  for 
integration  by  parts  which  will  affect  the  accuracy  in  determining  the 
instability.  It  has  been  shown  that  in  a  three-dimensional  case,  the  use 
of  tensorial  approach  is  more  efficient  than  the  vector  notations. 

The  viscosity  effect  is  small  as  expected  for  the  uniform  cylinder. 
However,  it  will  contribute  greatly  if  the  geometry  becomes  irregular. 

A  new  vortical  instability  integral  is  derived,  which  represents 
the  classical  hydrodynamic  instability.  This  avoids  solution  of  the 
three-dimensional  Orr-Summerfeld  equation.  A  simple  case  of  hyperbolic 
tangent  velocity  profile  is  used  for  shear  layer  instability  calculations. 
In  the  case  of  an  acoustic-coupled  vortical  instability,  it  is  noted  that 
no  contribution  results  from  tangential  modes. 

In  this  study,  we  utilize  an  analytical  form  of  classical  normal 
modes  for  the  cylindrical  acoustic  field.  For  irregular  geometeries,  how¬ 
ever,  we  must  resort  to  the  standard  eigenvalue  analysis  to  obtain  the 
acoustic  and  vortical  mode  shapes  and  frequencies.  With  such  data  we 
return  to  the  stability  integral  and  proceed  Identically  as  shown  in  the 
simple  examples  given  in  this  paper. 
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